Dislocation identification from orthogonal reconstructions¶
This notebook demonstrates a workflow to identify dislocations using orthogonal reconstructions. Only the reconstructed density and phase are required as inputs.
Overview:
- Load orthogonal reconstruction
- Extract amplitude (density) and phase
- Compute dislocation signature from phase discontinuities and density features
- Visualise candidate dislocation locations
Notes: Keep the input reconstructions aligned and sampled consistently across axes.
Load dataset¶
In [1]:
1
Out[1]:
1
In [2]:
import numpy as np
from matplotlib import rcParams
import cdiutils
from cdiutils.analysis import dislo
from cdiutils.utils import fill_up_support, nan_to_zero, zero_to_nan
from cdiutils.io.vtk import save_as_vti
# Plotting configurations
rcParams.update(
{
# #"font.family": "Liberation Serif", # Use for titles, labels, and body text
# # "font.size": 14, # General font size
# # "axes.labelsize": 14, # X and Y labels font size
# # "xtick.labelsize": 14, # X-axis tick labels
# # "ytick.labelsize": 14, # Y-axis tick labels
# # "legend.fontsize": 14, # Legend font size
# # "axes.titlesize": 14, # Title font size
"savefig.bbox": "tight",
# 'axes.linewidth': 2.0, # axis box (spines)
# 'grid.linewidth': 2.0,
# 'grid.linestyle': 'dotted',
"grid.alpha": 0.5,
"xtick.major.width": 0.5,
"ytick.major.width": 0.5,
"xtick.major.size": 0.5,
"ytick.major.size": 0.5,
"font.weight": "bold",
"axes.titleweight": "bold",
"axes.labelweight": "bold",
"figure.dpi": 150,
}
)
In [ ]:
In [3]:
# path_to_data = "path/to/your/dataset.npy"
path_to_data = "C:\\Users\\aboud\\data\\S457_amp-phase_bulk.npz"
data = np.load(path_to_data)
density = data["amp"]
phase = -data[
"phase"
] # the minus since this is the raw phase of pynx reconstruction
voxel_size = [float(i) for i in data["voxel_size"]] # in nm
new_mask = nan_to_zero(density) / np.nanmax(density) > 0.05
new_mask = fill_up_support(fill_up_support(new_mask))
density *= new_mask
phase *= new_mask
obj = density * np.exp(1j * phase)
In [4]:
figure, axes = cdiutils.plot.plot_volume_slices(
zero_to_nan(phase),
plot_type="contourf",
title="Phase field",
figsize=(10, 4),
label_size=12,
cmap="jet",
vmin=-np.pi,
vmax=np.pi,
)
# figure.savefig("original_phase.png")
In [5]:
from cdiutils.interactive import Plotter
Plotter((density), plot="3D")
Error updating plot: No surface found at the given iso value. Error updating plot: No surface found at the given iso value.
Out[5]:
<cdiutils.interactive.plotter.Plotter at 0x230f2b19970>
In [6]:
1
Out[6]:
1
compute the minimum gradient of the phase field¶
In [7]:
support_strain, strain_min = dislo.map_min_gradient(
obj=obj,
voxel_size=voxel_size,
save_filename_vti="dislo_fig/ortho_gradient_map.vti",
plot_debug=False,
path_to_save="",
save_plot=False,
)
# apply a threshoold
stain_threshold = 0.35
strain_min[strain_min < 0.3] = 0
Analysis¶
cleaning¶
here you can cut the segmented phase jump as there are noise espacilly near the surface and the orthogonalisation procesure.
The best is to compute this parametre before orthogonalisation then orthogonolise it. It is much better and cleaner
In [8]:
# use this plot and the vti to see if the clustering later works or no and adjust the parameters.
# figure, axes=cdiutils.plot.plot_volume_slices(
# zero_to_nan(phase)[:,:,40:],plot_type='contourf',
# title="Phase field", figsize=(10,4),label_size=12,
# cmap='jet',
# vmin=-np.pi, vmax=np.pi
# );
In [9]:
# if the strain_min is computed before the orthogonalisation then orthogonalised itself, much better and cleaner (which is my case)
strain_min = data["strain_amp"]
In [10]:
# strain_min[:35]=0
# strain_min[:,:35]=0
# strain_min[:,:,:40]=0
# strain_min[60:]=0
# strain_min[:,55:]=0
# strain_min[:,:,55:]=0
In [11]:
from cdiutils.interactive.multiviewer_3d import MultiVolumeViewer
In [12]:
dict_data = {
"density": density,
"phase": phase,
"dislo": strain_min,
}
viewer = MultiVolumeViewer(
dict_data,
figsize=(7, 9),render_workers=4,render_in_flight=4,rendering_mode='fast'
)
viewer.show()
clusering the voxels that have a phase jump¶
In [13]:
final_labeled_clusters, num_final_clusters = dislo.clusters_dislo_strain_map(
data=strain_min,
amp=density,
phase=phase,
save_path="step1",
voxel_sizes=tuple(voxel_size),
threshold=stain_threshold,
min_cluster_size=2,
distance_threshold=2.0,
cylinder_radius=1,
num_spline_points=5000,
smoothing_param=4,
eps=2.5,
min_samples=10,
save_output=True,
debug_plot=False,
)
# Save or visualize the circular mask and polar angles#
dict_to_vti = {
"density":nan_to_zero(density),
"phase":nan_to_zero(phase),
"dislo":final_labeled_clusters,
"strain_min":strain_min,
}
save_as_vti(output_path="dislo_fig/Step1.vti",voxel_size=tuple(voxel_size),**dict_to_vti)
Number of clusters identified: 3 Filtered clusters: [3] Cylindrical mask constructed.
unify the clusters related to dislocation line and create a cylinder based on the cluster¶
In [14]:
threshold_dislo_cluster = 0.5
unique_labels = np.unique(final_labeled_clusters)
print(unique_labels)
points = dislo.extract_structure(
final_labeled_clusters, threshold_dislo_cluster
)
# Fit a line to the extracted points
centroid, direction = dislo.fit_line_3d(points)
# Generate a filled cylinder in the same shape as the original volume
filled_cylinder_volume = dislo.generate_filled_cylinder(
shape=final_labeled_clusters.shape,
centroid=centroid,
direction=direction,
radius=1,
height=100,
step=1,
)
selected_dislocation_data = filled_cylinder_volume * new_mask
# Save or visualize the circular mask and polar angles#
dict_to_vti = {
"density":nan_to_zero(density),
"phase":nan_to_zero(phase),
"dislo":selected_dislocation_data,
"strain_min":strain_min,
}
save_as_vti(output_path="dislo_fig/Step2.vti",voxel_size=tuple(voxel_size),**dict_to_vti)
# np.savez_compressed("step1_2_dislo",**dict_to_vti)
[0 1]
In [15]:
dict_data = {
"density": density,
"phase": phase,
"strain_min": strain_min,
"dislo": strain_min,
}
viewer = MultiVolumeViewer(
dict_data,
figsize=(7, 9),render_workers=4,render_in_flight=4,rendering_mode='fast'
)
viewer.show()
plot phase around dislocation¶
In [16]:
radius_1 = 7 # 1 -> 15 nm
dz = 2
dr = 2
selected_point_index = 0
threshold_dislo_cluster = 0.5
points = dislo.extract_structure(
selected_dislocation_data, threshold_dislo_cluster
)
centroid, exp_direction = dislo.fit_line_3d(points)
# the direction is the direction of the dislocation line in the current frame after segmentation
In [17]:
# the displacement_vectors contain the vector of the voxel in the ring choosed and the centre of the dislocation in the current ref frame
(
phase_ring_3d,
angle_ring_3d,
circular_mask,
displacement_vectors,
direction,
) = dislo.plot_phase_around_dislo(
density,
phase,
selected_dislocation_data,
radius_1,
dr,
centroid,
exp_direction,
slice_thickness=dz,
voxel_sizes=voxel_size,
selected_point_index=selected_point_index,
save_vti=True,
fig_title=None,
save_path="dislo_fig/step3",
)
checking the results directly in notebook using MultiVolumeViewer
In [18]:
dict_data ={
"density": density,
"phase": phase,
"strain_min": strain_min,
"dislo": selected_dislocation_data,
"phase_ring_3d": phase_ring_3d,
"angle_ring_3d": angle_ring_3d,
"circular_mask": circular_mask,
}
viewer = MultiVolumeViewer(
dict_data,
figsize=(7, 12),render_workers=4,render_in_flight=4,rendering_mode='fast'
)
viewer.show()
In [19]:
phase_ring_3d = nan_to_zero(phase_ring_3d)
In [20]:
results___ = dislo.dislo_process_phase_ring(
angle_ring_3d,
phase_ring_3d,
displacement_vectors,
factor_phase=+1.0,
plot_debug=True,
save_path="dislo_fig/step4.pdf",
)
Raw angle : -179.34344 179.24577 Raw phase : -172.47305665095348 173.66661009827837 unwrapped phase : -190.2186981806769 312.03755268956263
CLean data in the jump part¶
In [21]:
conditions = ( # ( (phase_ring_3d-angle_ring_3d < 2.6) & (phase_ring_3d-angle_ring_3d >= -3.4) ) &
(angle_ring_3d < 0.52) & (angle_ring_3d > 0.3)
)
phase_ring_3d_clean = np.array(phase_ring_3d)
angle_ring_3d_clean = np.array(angle_ring_3d)
phase_ring_3d_clean[conditions] = 0
angle_ring_3d_clean[conditions] = 0
In [22]:
# (angle_raw, phase_raw, angle_final, phase_final, phase_smooth,phase_sinu, displacement_vectors_ring_sorted, displacement_vectors_final)
(_, _, angle_final, phase_final, _, _, _, displacement_vectors_final) = (
dislo.dislo_process_phase_ring(
angle_ring_3d_clean,
phase_ring_3d_clean,
displacement_vectors,
factor_phase=-1.0,
plot_debug=True,
save_path="dislo_fig/step4beta.pdf",
)
)
Raw angle : -179.34344 179.24577 Raw phase : -170.63139793033577 173.33815041665468 unwrapped phase : -51.11729422932207 312.03755268956263
Identification of burger vector from theoritical model¶
General parametre definition¶
In [23]:
# must be known for the current orthogonalised frame
Rotation_matix_cristallographic_to_crystal_final = np.array(
[
[-0.75403343, 0.64962679, 0.09704953],
[-0.3200682, -0.49242258, 0.80936787],
[0.57357644, 0.57922797, 0.57922797],
]
)
Rotation_matix_crystal_to_cristallographic_final = np.linalg.inv(
Rotation_matix_cristallographic_to_crystal_final
)
In [24]:
# the displacement_vectors contain the vector of the voxel in the ring choosed and the centre of the dislocation in the current ref frame
vx_example = displacement_vectors_final[..., 0]
vy_example = displacement_vectors_final[..., 1]
vz_example = displacement_vectors_final[..., 2]
## Transform vector to crystallographic basis using R_final
vx_cryst, vy_cryst, vz_cryst = (
dislo.transform_known_vector_to_crystallographic(
vx_example,
vy_example,
vz_example,
Rotation_matix_crystal_to_cristallographic_final,
)
)
vx_cryst, vy_cryst, vz_cryst = dislo.normalize_vectors_3d(
vx_cryst, vy_cryst, vz_cryst
)
vx_cryst, vy_cryst, vz_cryst = (
vx_cryst / vx_cryst,
vy_cryst / vx_cryst,
vz_cryst / vx_cryst,
)
exp_vect_Crystallographic = np.stack([vx_cryst, vy_cryst, vz_cryst]).T
In [25]:
exp_direction_crystallo = (
exp_direction @ Rotation_matix_cristallographic_to_crystal_final
)
# take the reference angle at polar angle theta 0
angle_closetozeo, index_0 = dislo.closest_to_zero_in_array(
angle_final,
)
# the vector of the voxel where the angle is close to 0
exp_vectdislo = exp_vect_Crystallographic[index_0]
# separate the parrallel (x) and perp (y) component of the vector compared to the dislocation direction exp_direction_crystallo
exp_xdislo = dislo.normalize_vector(
dislo.project_vector(exp_vectdislo, exp_direction_crystallo)
)
exp_ydislo = dislo.normalize_vector(
np.cross(exp_xdislo, exp_direction_crystallo)
)
angle_final *= np.pi / 180
phase_final *= np.pi / 180
In [26]:
a_0_pt = 3.924 / 10 # lattice paramettre for platinum in nm
d_0pt = a_0_pt / 3**0.5 # d-spacing for bragg 111
t = np.array(
[8.0, 3.0, 8.0]
) # must be close to exp_direction_crystallo (direction of the dislocation in crystallographic frame)
G_exp = np.array([1.0, -1.0, 1.0]) * (2 * np.pi / a_0_pt)
comparison to theo¶
In [27]:
# compute the theoritical prediction of the dislocation in infinite medium aproximation for the 3 possible burger vector
b_cases = np.array(
[
[-0.196, 0.196, 0.0],
[-0.196, 0.0, -0.196],
[0.0, 0.196, -0.196], # (-b)
# [ 0.196, -0.196, 0. ], [ 0.196, 0. , 0.196], [ 0. , -0.196, 0.196], #(+b)
]
)
b_cases_str = np.array(
[
str((i.astype(int))).replace("[", "").replace("]", "").replace(" ", "")
for i in np.round(b_cases / (0.5 * a_0_pt))
]
)
phase_theo_3_cases = [] # phase for the three possible fcc burger vector
for i_b in range(len(b_cases)):
print(f" ++++++++++++++++ b: {b_cases[i_b]} +++++++++++++++++++++")
predicted_phase = dislo.dislo_phase_model(
angle_final,
t,
G_exp,
b_cases[i_b],
align_theta=exp_xdislo,
print_debug=True,
)
predicted_phase, pred_fit, _, f_linear_sim, coeffs_linear = (
dislo.decompose_experimental_phase(angle_final, predicted_phase)
)
phase_theo_3_cases.append(predicted_phase)
++++++++++++++++ b: [-0.196 0.196 0. ] +++++++++++++++++++++ the bper is off by 139.7797395336527 ° from the experimental reference b_perp: [-0.13877372 0.21745985 0.05722628], b_perp_norm: 0.26423811424905724 b_screw: [-0.05722628 -0.02145985 -0.05722628] b_screw_norm: -0.08372705045624257 Rotation matrix R: [[-0.52518435 0.82296929 0.21657087] [-0.5069794 -0.5069794 0.69709668] [ 0.68348613 0.2563073 0.68348613]] G_theo: [-18.11912406 11.16204788 17.78418381] ++++++++++++++++ b: [-0.196 0. -0.196] +++++++++++++++++++++ the bper is off by 108.15375256851857 ° from the experimental reference b_perp: [-0.01287591 0.06867153 -0.01287591], b_perp_norm: 0.07104475817562882 b_screw: [-0.18312409 -0.06867153 -0.18312409] b_screw_norm: -0.26792656145997623 Rotation matrix R: [[-0.18123663 0.96659535 -0.18123663] [-0.70710678 0. 0.70710678] [ 0.68348613 0.2563073 0.68348613]] G_theo: [-21.28130564 0. 17.78418381] ++++++++++++++++ b: [ 0. 0.196 -0.196] +++++++++++++++++++++ the bper is off by 76.52487147248954 ° from the experimental reference b_perp: [ 0.05722628 0.21745985 -0.13877372], b_perp_norm: 0.26423811424905724 b_screw: [-0.05722628 -0.02145985 -0.05722628] b_screw_norm: -0.08372705045624257 Rotation matrix R: [[ 0.21657087 0.82296929 -0.52518435] [-0.69709668 0.5069794 0.5069794 ] [ 0.68348613 0.2563073 0.68348613]] G_theo: [-18.11912406 -11.16204788 17.78418381]
In [28]:
dislo.plot_phase_data_comparison_exp_to_theo(
angle_final,
phase_final,
phase_theo_3_cases,
b_cases_str,
save_path=None,
filter_low_freq=False,
marker_size=5,
figsize=(9, 5),
font_size=12,
show_slope=True,
ncol=4,
fix_exp_slope=0.98,
)
comparison to simulation ....¶
In [ ]: